; +
; NAME:
; gauge_flux_mozart_3D_regrid
;
; PURPOSE:
; Reformats it to mozart grid and writes the file
;
; CATEGORY:
; regridding
;
; CALLING SEQUENCE:
; GAUGE_FLUX_MOZART_3D_REGRID, q_lon, q_lat, q, date, flux_ncname, infilename
;
; INPUTS:
; q_lon = Input longitudes
; q_lat = Input latitudes
; q = Input flux
; date = date from input flux file
; datesec = datesec from input flux file
; flux_ncname = Flux name output from gauge_flux_extract_file_info
; infilename = Input filename
; varinfo = variable information from gauge_flux_extract_file_info
; species = species as string
;
; OUTPUTS:
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
; Regrids horizontally each level in a loop onto MOZART grid
; Regrids vertically onto the Pressure levels of a MOZART MERRA file
; Does the above for a single year and then replicates it for 10 years starting in 2000
;
; MODIFICATION HISTORY:
;   Written by:  Ann Stavert 2014
;   Modified from gauge_flux_mozart_2D_regrid
;
; -

Pro GAUGE_FLUX_MOZART_3D_REGRID, q_lon, q_lat, q, date, datesec, flux_ncname, infilename, varinfo, species, $
  TITLE = Title, PS_FILE = PS_File
  
  ; need to extract the level variable
  gauge_flux_extract_file_info, INFILENAME, DimNames, DimSizes,  VarInfo, PS_NCNAME = PS_ncname, GLOBAL_ATTS= global_atts
  
  if keyword_set(PS_File) then q_PS=mr_ncget('PSURF', filename=PS_FILE) else q_PS=mr_ncget(PS_ncname, filename=InFileName)
  Varinfo_PS = varinfo[strjoin(PS_ncname)]
  PS_attnames = Varinfo_PS.AttName
  PS_atts = VarInfo_PS.Attributes
  PS_units = PS_atts[where(strmatch(PS_attnames, 'units', /FOLD_CASE))]
  
  
  ;Create some info for the file
  If keyword_set(TITLE) eq 0 then Title = Species + '_emissions'
  Author='Ann Stavert'
  
  ;Get flux data
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  lon_ncname='lon'
  lat_ncname='lat'
  
  
  ; Extract the data from the flux file
  var_hyam =mr_ncget('A', filename=infilename) ;var_hyam =mr_ncget('hyam', filename=infilename)
  var_hybm =mr_ncget('B', filename=infilename) ;var_hybm =mr_ncget('hybm', filename=infilename)
  
  
  
  ; Orientate gridded data, so that latitudes are ordered from -90 to 90, and
  ;   longitudes are ordered from 0 to 360.
  lattemp=q_lat
  lontemp=q_lon
  mr_grid_orient, q_lon, q_lat, q
  mr_grid_orient, lontemp, lattemp, q_PS
  
  ;change date ('MM') to Julian date
  time=julday(fix(date), 15, 2000,0,0,0)
  
  ; generate pressure field for input grid
  var_P=fltarr(n_elements(q_lon), n_elements(q_lat), n_elements(var_hyam), 12)
  for ti=0, 11 do var_P[*, *, *, ti]=mr_hybrid_coords(var_hyam, var_hybm, q_PS[*, *, ti])
  
  
  ;Get MOZART grid
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  MOZART_fname= '/data/shared/GAUGE/MOZART/Mean_MERRA_Monthly_PS.nc'
  MOZART_T_fname= '/data/shared/GAUGE/MOZART/Mean_MERRA_Monthly_T.nc'
  lon=mr_ncget('lon', filename=MOZART_fname)
  lat=mr_ncget('lat', filename=MOZART_fname)
  lev=(mr_ncget('lev', filename=MOZART_fname))[16:-1]
  hyam=(mr_ncget('hyam', filename=MOZART_fname))[16:-1] ;Only keep lowest 54 levels
  hybm=(mr_ncget('hybm', filename=MOZART_fname))[16:-1] ;Only keep lowest 54 levels
  PS=mr_ncget('PS', filename=MOZART_fname)
  T=(mr_ncget('T', filename=MOZART_T_fname))[*, *, 16:-1, *]
  
  PSa=fltarr(n_elements(lon), n_elements(lat), 12)
  Ta=fltarr(n_elements(lon), n_elements(lat), n_elements(hyam), 12)
  nYears=n_elements(PS[0, 0, *])/12
  for mi=0, 11 do begin
    PSa[*, *, mi]+=total(PS[*, *, indgen(nYears)*12 + mi], 3)
   ; maxmin, PS[*, *, indgen(nYears)*12 + mi]
    Ta[*, *, *, mi]+=total(T[*, *, *, indgen(nYears)*12 + mi], 4)
    ; maxmin, T[*, *, indgen(nYears)*12 + mi]
  endfor
  PS=temporary(PSa)/float(nYears)
  T=temporary(Ta)/float(nYears)
  
  ;Define output array
  var_regrid_horizontal=fltarr(n_elements(lon), n_elements(lat),n_elements(var_hybm), 12)
  var_p_regrid_horizontal=fltarr(n_elements(lon), n_elements(lat),n_elements(var_hybm), 12)
  var_regrid=fltarr(n_elements(lon), n_elements(lat),n_elements(hybm), 12)
  
  ; generate pressure field for mozart grid
  P=fltarr(n_elements(lon), n_elements(lat), n_elements(hybm), 12)
  for mi=0, 11 do begin
    P[*, *, *, mi]=mr_hybrid_coords(hyam, hybm, PS[*, *, mi], P0=100000.)
  endfor
  
  lonII=mr_interp_indices(q_lon, lon)
  latII=mr_interp_indices(q_lat, lat)
  
  
  ; Regrid 2D slices of the variable of interest by iterating through each pressure level and time
  for i = 0L, n_elements(var_hybm)-1 do begin ; loop through each pressure level
  
    for ti=0, n_elements(time)-1 do begin ; loop through each time
      ;Regrid
      var_regrid_horizontal[*, *, i, ti]=bilinear(q[*, *, i, ti], lonII, latII) ; regrid the variable of interest
      var_p_regrid_horizontal[*, *, i, ti]=bilinear(var_p[*, *, i, ti], lonII, latII) ; regrid the pressure
      print, 'Horizontal regridding ' + string((float(ti) + n_elements(time)*i )/float(n_elements(time)*n_elements(var_hybm))*100., format='(I)'), '% complete. ' + $
        string(round(mean(var_regrid_horizontal[*, *,i, ti])/mean(q[*, *,i, ti])-1.)*100., format='(I)'), '% error'
        
    endfor ;end ti loop
    
  endfor ;end i loop
  
  
  ; Now have the 3D field on the lat/lon of the MERRA
  ; Now need to interpolate each column from the pressures of the intital grid to that of the MERRA grid
  ; Use the MERRA_P and var_P variables
  
  for lonI=0, n_elements(lon)-1 do begin
    for latI=0, n_elements(lat)-1 do begin
      for mi=0, 11 do begin
        var_regrid[lonI, latI, *, mi]=interpol(var_regrid_horizontal[lonI, latI, *, mi], $
          var_p_regrid_horizontal[lonI, latI, *, mi], P[lonI, latI, *, mi])
      endfor
    endfor
  endfor
  
  
  ;convert units to VMR
  Varinfo_q = varinfo[flux_ncname]
  q_attnames = Varinfo_q.AttName
  q_atts = VarInfo_q.Attributes
  q_units = q_atts[where(strmatch(q_attnames, 'units', /FOLD_CASE))]
  
  
  case q_units of
    "CM-3": var_regrid=var_regrid*1.e6*1.381e-23*T/P ; NB: uses boltzmanns constant = gas const (R) /avagadros number
    "cm-3": var_regrid=var_regrid*1.e6*1.381e-23*T/P
    "cm^-3": var_regrid=var_regrid*1.e6*1.381e-23*T/P
    "s^-1": var_regrid=var_regrid*1.e6*1.381e-23*T/P
    "cm-3 s-1": var_regrid=var_regrid*1.e6*1.381e-23*T/P
    else: begin
    
      Convert_To_MoleculesCm2Sec, var_regrid, q_out, q_Units
      var_regrid = q_out
      
    end
  endcase
  
  
  ;Write out the file for a number of years
  ; Make a new date stamp for each year
  startY=1988
  EndY=2012
  nYears=EndY - StartY +1
  
  outputfiles = !null
  
  
  for i = 0L, nYears-1 do begin
    time_out=timegen(12, start=julday(1, 15, startY+i, 0), unit='monthly')
    
    If total(strmatch(['CO2','N2O'], Species, /FOLD_CASE)) ne 1 then begin
      if keyword_set(NCDirOut) eq 0 then ncDirOut = DIALOG_PICKFILE(Path = '/data/shared/GAUGE/MOZART/', TITLE= 'Select output directory', /DIRECTORY)
    endif else  ncDirOut = '/data/shared/GAUGE/MOZART/' + SPECIES +'/'
    
    cd, ncDirOut, current = current
    if file_test( strtrim(StartY + i, 2), /DIRECTORY) eq 0 then File_MKDIR, strtrim(StartY + i, 2)
    cd, current
    
    
    ;ncName = strjoin((strsplit(file_basename(INfilename, '.nc'), "_", /EXTRACT)),'.') + '.' + strtrim(startY+i, 2) + '.nc'
    ncname = SPECIES + '0001.nc'
    
    ;Now write out MOZART netCDF file
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    fID=ncdf_create(ncDirOut + strtrim(StartY + i, 2) +'/' + ncName, /clobber)
    
    TimeDID  = NCDF_DIMDEF(fID, 'time', /unlimited)
    LonDID  = NCDF_DIMDEF(fID, 'lon', long(n_elements(lon)))
    LatDID  = NCDF_DIMDEF(fID, 'lat', long(n_elements(lat)))
    LevDID  = NCDF_DIMDEF(fID, 'lev', long(n_elements(lev)))
    
    ; Create a variable to hold the data
    lonID = NCDF_VARDEF(fId,'lon', [LonDID], /DOUBLE)
    latID = NCDF_VARDEF(fId,'lat', [LatDID], /DOUBLE)
    levID = NCDF_VARDEF(fId,'lev', [LevDID], /DOUBLE)
    hyamID = NCDF_VARDEF(fId,'hyam', [LevDID], /DOUBLE)
    hybmID = NCDF_VARDEF(fId,'hybm', [LevDID], /DOUBLE)
    dateID = NCDF_VARDEF(fId,'date', [TimeDID], /LONG)
    timeID = NCDF_VARDEF(fId,'time', [TimeDID], /DOUBLE)
    varID = NCDF_VARDEF(fId, flux_ncname, [LonDID, LatDID, levDID, TimeDID], /DOUBLE)
    psID = NCDF_VARDEF(fId, 'PS', [LonDID, LatDID, TimeDID], /DOUBLE)
    
    ;Create some global attributes
    NCDF_ATTPUT, fid, /Global, "Title", strjoin(title)
    NCDF_ATTPUT, fid, /Global, "Author", author
    NCDF_ATTPUT, fid, /Global, "Created", systime()
    
    ; Create some attributes (about our variables)
    NCDF_ATTPUT, fId, lonId, "long_name", "longitude"
    NCDF_ATTPUT, fId, lonId, "units", "degrees_east"
    
    NCDF_ATTPUT, fId, latId, "long_name", "latitude"
    NCDF_ATTPUT, fId, latId, "units", "degrees_north"
    
    NCDF_ATTPUT, fId, levId, "long_name", "level"
    NCDF_ATTPUT, fId, levId, "units", "hPa"
    
    NCDF_ATTPUT, fId, hyamId, "long_name", "Hybrid A coordinate"
    NCDF_ATTPUT, fId, hyamId, "units", "Pa"
    
    NCDF_ATTPUT, fId, hybmId, "long_name", "Hybrid B coordinate"
    NCDF_ATTPUT, fId, hybmId, "units", " "
    
    NCDF_ATTPUT, fId, varId, "long_name", flux_ncname
    NCDF_ATTPUT, fId, varId, "units", "VMR"
    
    NCDF_ATTPUT, fId, psId, "long_name", "Surface Pressure"
    NCDF_ATTPUT, fId, psId, "units", "Pa"
    
    NCDF_ATTPUT, fId, TimeID, "long_name", "Time"
    NCDF_ATTPUT, fId, TimeID, "units", "days since 0000-01-01 00:00:00"
    
    NCDF_ATTPUT, fId, DateID, "long_name", "Date"
    NCDF_ATTPUT, fId, DateID, "units", "YYYYMMDD"
    
    ; Leave definition mode and enter data write mode
    NCDF_CONTROL, fId, /ENDEF
    
    ; Write the data
    NCDF_VARPUT, fId, lonId, lon
    NCDF_VARPUT, fId, latId, lat
    NCDF_VARPUT, fId, levId, lev
    NCDF_VARPUT, fId, hyamId, hyam
    NCDF_VARPUT, fId, hybmId, hybm
    
    NCDF_VARPUT, fId, dateID, mr_mozdate(time_out)
    NCDF_VARPUT, fId, timeID, time_out - julday(1, 121, 1, 0)
    NCDF_VARPUT, fId, varId,  var_regrid
    NCDF_VARPUT, fId, psId,  PS ; Surface P
    
    ; Done
    NCDF_CLOSE, fId
    
    print, 'Wrote file : ' + ncdirout + strtrim(StartY + i, 2) + '/' + ncName
    
    Outputfiles = [OutputFiles, ncdirout + strtrim(StartY + i, 2) + '/' + ncName]
  endfor
  
  
  ; Do appending now!
  ; Now go back and append a data point at the start and end of each file
  AppendStartAndEnd_Mozart_NC, OutputFiles, SPECIES=SPECIES, VARIABLE = flux_ncname, UNITS = 'VMR', TITLE = Title
  
  
END